A  Variable-metric  Variant  of  the 
Karmarkar  Algorithm  for  Linear  Programming 


J.E.  Dennis,  Jr.1 

A.M.  Morshedi2 
and 

Kathryn  Turner3 

Technical  Report  86-13,  June  1986 
Revised  July  1986,  September  1986,  and  January  1987 


'Mathematical  Sciences  Department,  Rice  University,  Houston,  Texas  77251-1892.  Research  sponsored  by  DOE  DE- 
AS05-82ER 13016  ARO  DAAG-29-83-K-0035,  AFOSR  85-0213 

^ hell  Development  Company,  Westhollow  Research  Center,  Houston,  Texas  77001 

Mathematical  Sciences  Department,  Rice  University,  Houston,  Texas  77251-1892.  Research  sponsored  by  ARO 
DAAG-29-83-K-0035,  AFOSR  85-0243,  Shell  Development  Company. 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

JAN  1987  2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-1987  to  00-00-1987 

4.  TITLE  AND  SUBTITLE 

A  Vairable-metric  Variant  of  the  Karmarkar  Algorithm  for  Linear 
Programming 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROIECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Computational  and  Applied  Mathematics  Department  ,Rice 

University, 6100  Main  Street  MS  134, Houston, TX, 77005-1892 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBIECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

18.  NUMBER  19a.  NAME  OF 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE 

unclassified  unclassified  unclassified 

34 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


A  Variable- metric  Variant  of  the  Karmarkar  Algorithm 

for  Linear  Programming 

J.E.  Dennis,  Jr.1,  A.M.  Morshedi2,  and  Kathryn  Turner3 

Technical  Report  86-13,  June  1986 
Revised  July  1986,  September  1986,  and  January  1987 


Abstract 


The  most  time-consuming  part  of  the  Karmarkar  algorithm  for  linear  pro¬ 
gramming  is  the  projection  of  a  vector  onto  the  nullspace  of  a  matrix  that 
changes  at  each  iteration.  We  present  a  variant  of  the  Karmarkar  algorithm  that 
uses  standard  variable-metric  techniques  in  an  innovative  way  to  approximate 
this  projection.  In  limited  tests,  this  modification  greatly  reduces  the  number  of 
matrix  factorizations  needed  for  the  solution  of  linear  programming  problems. 
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1.  Introduction:  the  K&rmark&r  algorithm 


Since  the  introduction  of  the  polynomial  time  algorithm  for  linear  program¬ 
ming  by  N.  Karmarkar  [1984],  there  has  been  considerable  interest  in  its  use,  and 
in  developing  strategies  for  modifying  the  algorithm  to  achieve  greater  computa¬ 
tional  efficiency.  We  present  a  variable-metric  variant  of  the  Karmarkar  algo¬ 
rithm.  The  algorithm  as  proposed  by  Karmarkar  employs  a  projective  transfor¬ 
mation  to  map  the  feasible  region  of  the  linear  program  onto  the  intersection  of  a 
simplex  with  an  affine  space,  and  then  solves  the  problem 

minimize  cTx 

subject  to  Ax  =  0 
T 

ex  =  n 
x  >  0, 

where  c,x,e  6  IRn,  A  €  IR'”*",  and  e=(l,  •  •  •  ,l)r  We  follow  the  suggestion 
of  Shanno  and  Marsten  [1985]  in  using  eTx  —  n,  rather  than  eTx  —  1,  as  used 
by  Karmarkar  [1984].  We  transform  a  standard  linear  programming  problem  to 
the  above  form  by  commonly  used  strategies  that  do  not  involve  a  projective 
transformation.  For  completeness,  we  include  these  strategies  in  Appendix  I.  We 
shall  also  assume  that  the  minimum  value  of  the  objective  function  cT x  is  zero, 
and  that  we  have  an  initial  feasible  point  x0>0.  A  means  of  getting  such  an  ini¬ 
tial  feasible  point  and  of  transforming  a  given  objective  function  to  an  objective 
function  whose  minimum  value  is  zero  is  also  included  in  Appendix  I. 

At  each  step  of  the  Karmarkar  algorithm,  the  current  iterate  (x1?  •  •  •  ,xn)T 
is  mapped  to  the  center  (1,  •  ■  ■  ,l)r  of  the  simplex  by  the  projective  transforma¬ 
tion 


x  — 


TD~\ 


■ D -1 


x, 


where  D  =  diag  (ij,  •  •  •  ,xn),  and  we  use  x  to  denote  the  transformed  variable. 
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A  step  is  taken  in  the  transformed  space  in  the  direction  of  the  negative  pro¬ 
jected  gradient  of  a  linear  function,  and  the  new  point  in  the  untransformed 
space  is  found  by  applying  the  inverse  of  the  above  projective  transformation 
(restricted  to  {i  €  Kn  :  x  >  0,  eTx  =  n}), 


The  Karmarkar  algorithm  is  : 


given  an  initial  feasible  point  x0  >  0 

begin 

x  z0 

while  c  T  x  is  too  large 

do 

D  =  diag  (®j,  •  •  •  ,z„), 


project  Dc  onto  the  nullspace  of  B  : 
Cp  =  [I-BT{BBT)-lB]Dc 


x  —  t  —  ac. 


x  — 


n 


t  T  Dx 


■Dx 


end  do 


end 


In  section  3  we  discuss  how  the  steplength  parameter  a  is  chosen  at  each  step. 
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2.  A  variable- metric  variant  of  the  Karmarkar  algorithm 


The  most  time-consuming  part  of  the  Karmarkar  algorithm  is  the  projection 
of  the  vector  Dc  onto  the  nullspace  of  a  matrix  B  that  changes  at  each  iteration. 
We  propose  to  approximate  the  direction  so  obtained,  using  secant  updates  to  the 
initial  matrix.  The  updating  strategies  we  have  adapted  for  use  in  this  context 
have  been  used  very  successfully  in  the  solution  of  nonlinear  equations  and  in  the 
minimization  of  nonlinear  functions.  By  employing  updates  we  reduce  the 
amount  of  work  required  at  each  step,  since  we  bypass  the  need  to  obtain  a 
matrix  factorization  at  every  step.  We  shall  approximate  the  direction 
Cp  =  [I-Bt{BBt)~1B]Dc  by 

cp  =  D-1D[I-BT{BBTYlB]DTc,  (2.1) 

where 


AD 

cTD~lD 


BDlD, 


and  D  is  a  nonsingular  approximation  to  D.  The  new  iterate  in  the  transformed 
space  is  given  by 


x  =  e  -aCp  . 

Before  addressing  the  way  we  obtain  the  matrix  D ,  we  note  that  any  direction  cp 
computed  in  this  manner  is  a  feasible  direction  for  the  linear  program  and  the 
negative  of  it  is  a  descent  direction  for  Karmarkar’s  potential  function, 

/(*)=  £log-^.  (2.2) 

Feasibility  of  the  direction  cp  follows  from  the  observation  that  cp  is  in  the 
nullspace  of  B.  Karmarkar’s  proof  of  the  convergence  of  his  algorithm  and  his 
polynomial  time  bound  are  based  on  achieving  at  least  a  constant  amount  of 
reduction  in  the  potential  function  (2.2)  at  each  step.  This  potential  function  is 
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transformable  by  the  projective  transformation  to  a  function  having  the  same 
form.  When  /  is  expressed  in  terms  of  the  transformed  variable,  we  have 


/(*)  =  £log 

1  =  1 


c  T  Pi 
i  i 


+  constant, 


so  that  in  the  transformed  space,  we  may  consider 

”  cT  Dx 

/(*)  =  £ 

!  =  1  X  i 

Now 


(2.3) 


V-  /  (i  )  —  — —  Dc  -  D  1e,  where  D  —  diag  (i  v  •  •  •  ,x  B)  , 
c1  Dx 

and 

Vs/  (e)  —  "  Dc  -  e, 

c  xc 

where  xc  denotes  the  current  iterate  in  the  untransformed  space.  The  direction 
-cp  is  a  non-ascent  direction  for  the  potential  function  in  the  transformed  space 
since 

Vf/  (c)r(-Cp)  =  j-^“  Dc  +  eJ 

=  — £—cTD[I-BT{BBTy1B]DTc 
c  xc 

<  0  . 


In  order  to  show  that  the  direction  -cp  is  a  descent  direction  for  the  poten¬ 
tial  function,  we  note  that  cp  7^  0  ,  since  there  exists  a  step  in  the  direction  -cp 
that  yields  a  reduction  in  the  potential  function,  and  that  -cp  can  fail  to  be  a 
descent  direction  only  if  D  c  is  orthogonal  to  the  nullspace  of  B.  But  if  D  c  is 
orthogonal  to  the  nullspace  of  B ,  then  there  exists  y  6  IRm  such  that 
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BTy  =  DT  c 

DT D~lBT y  =  DT c ,  and  since  £>  is  nonsingular, 

BTy  =  Dc, 

which  contradicts  the  fact  that  cp  7^  0  . 

Our  modification  of  the  Karmarkar  algorithm  requires  a  factorization  of  the 
initial  matrix  ( AD0)(AD0)T  only.  Thereafter,  as  is  the  usual  approach  in  compu¬ 
tations  for  large  problems,  update  vectors  are  stored,  and  the  projection  compu¬ 
tation  is  done  cheaply,  requiring  only  solutions  of  linear  systems  from  the  initially 
factored  matrix,  scalar  products,  and  sums  of  vectors. 

A 

We  begin  with  D0  —  D0,  a  diagonal  matrix  whose  diagonal  elements  are  the 
components  of  the  initial  feasible  point  x0.  We  may  obtain  rank-one  secant 
updates  to  this  matrix  through  the  use  of  one  of  two  nonlinear  functions  whose 
gradients  and  Hessians  involve  the  matrix  D.  The  first  and  simplest  of  these  is  a 
logarithmic  barrier  function,  and  the  second  is  the  potential  function  we  have 
already  defined.  Let  us  first  examine  the  logarithmic  barrier  function. 

Applying  the  barrier  transformation  to 

min  cT  x 

X 

subject  to  Ax  =  0  (24) 

eT  x  =  n 
x  >  0 


gives 


F(x)  =  c  T x  -  fi  log  x, 

»'=i 

subject  to  Ax  =  0 


min 

X 


e  x  —  n  , 


(2.5) 


the  inequality  constraints  having  been  replaced  by  a  term  in  the  objective  func¬ 
tion.  It  is  well  known  (see,  for  example,  Gill,  Murray,  Saunders,  Tomlin,  and 
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Wright,  [1985])  that  under  mild  hypotheses  the  solution  to  (2.5)  converges  to  the 
solution  of  (2.4)  as  /*— ►  0.  It  is  not  our  intention  to  solve  problem  (2.5)  for  any 
value  of  ji;  we  intend  merely  to  use  the  logarithmic  barrier  function,  whose  gra¬ 
dient  and  Hessian  involve  the  matrix  D  ,  in  order  to  obtain  an  approximation  to 
this  matrix.  We  have 

VxF(x)  =  c  -  nD~lt, 

where  D  —  diag  (x1?  •  ■  •  ,xn)  and  c  =  (1,  •  •  ,1)^;  and 

Vz2F(x)=/iP-2. 

Analogous  to  the  idea  of  a  secant  approximation  to  the  derivative  of  a  function  of 
one  variable  as  the  ratio  of  the  change  in  function  values  to  the  change  in  the 
independent  variable,  a  secant  approximation  to  a  Hessian  is  obtained  by  requir¬ 
ing  the  new  approximation  to  the  Hessian  acting  on  the  step  in  the  independent 
variable  to  produce  the  difference  in  gradient  vectors  at  the  new  and  previous 
points.  Using  the  subscript  +  to  refer  to  the  new  iterate  and  the  subscript  c  to 
refer  to  the  previous  (current)  iterate,  we  may  write  the  secant  equation  based  on 
the  logarithmic  barrier  function  as 

fiD+2{x+-xc)  =  (c  -  pD?e)-{c  -nD~le) 

D+2{x+-xc)  =  {D~l  -  D+1  )e. 

We  compute  the  righthand  side  exactly.  Letting  y b  denote  (D^1  -  D+1  )e  and, 
not  wishing  to  require  the  approximation  D  to  be  symmetric,  replacing  D  2  with 
D+D+,  we  require  the  approximation  to  satisfy 

{x+-xc)=D+blyb. 

An  approximation  based  on  the  potential  function  (2.2)  requires  a  very  simi¬ 
lar  secant  equation  to  be  satisfied.  Recall  that  the  potential  function  is 
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so  that 


/(*)  =  E1o8 


t 

C  1  X 


i= 1 


» 


v,/(*) 

v,V(x) 

and  the  secant  equation  is 


n 


T 

c  1  x 


c  -  D  le 


—  n-2 


/>-*- 


n 

(cri)2 


cc 


P+2  "  7  r"  ~r2  cc  r)  )  =  ~~T — c  -  D+le\  -  -j—c  -Dcle 

(c  1  x  ,  Y  CJI,  c  ar. 


Moving  the  second  term  on  the  lefthand  side  of  the  equation  to  the  righthand 
side  and  collecting  terms  gives 


Df  (*+  -x,)  =  {D-'-Dl1  )e 


n(cTzc-cTx+)2 

(cTx+f(cTxc) 


(2.7) 


Now  letting  yp  denote  the  righthand  side  of  the  above  equation,  and  again 

An  a  aw 

replacing  D+  with  D+D+ ,  we  require  the  approximation  to  satisfy 

(x+-xc)  =  b+i>lyv  ■ 


Recall  that  the  matrix  we  are  trying  to  approximate  is  a  known  diagonal 
matrix.  A  third  and  even  simpler  alternative  requirement  that  can  be  imposed  on 
the  approximation  is  that  it  satisfy 

D+2(x+-xc)  =  Df{x+-xc)  .  (2.8) 

We  shall  use  yt  to  denote  Z?+2(x+-zc)  and  express  this  requirement  as 

(i-f-x,.)  =  b+D+yt  . 


To  obtain  an  updating  formula  based  on  the  barrier  function  or  on  the 
potential  function  or  on  the  true  matrix  D+  ( y  denotes  either  yb  or  yp  or  $^,the 
righthand  side  of  equation  (2.6)  or  (2.7)  or  (2.8),  respectively),  we  follow  Dennis 
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and  Schnabel  [1983].  For  completeness,  we  include  the  development  of  the 
updating  formula  here.  We  express  the  requirement  ( x+-xc )  =  D+D+y  as 


(*+-**)  = 

v  =  £>ly 


(2.9) 


for  some  vector  v  6  1RB  •  Now  we  use  the  Broyden  updating  formula  to  obtain 
D+  as  the  closest  matrix  to  Dc  consistent  with  satisfying  the  first  part  of  (2.9) : 

( x+-xc-Dcv)vt 


D+  =  Dc  + 

The  second  part  of  (2.9)  requires  that 


T 

vx  v 


AT  AT  ,  ( x+-xe-Dev)Ty 
v  —  D+y  =  Dc  y  + - - - -  v 


T 

v  v 


and  this  can  be  satisfied  only  if 


«  -  *  bjv 


(2.10) 


(2.11) 


(2.12) 


for  some  k  €  K.  Now  substituting  (2.12)  into  (2.11)  and  simplifying  gives 

{x+-xc)Ty 


k2  = 


yTDcDjy 


Assuming  {x+-xc)T y  >  0,  we  choose  the  positive  root,  so  that 

j  (x+-xc)Ty  ] 

v  = 


(bJy)T(D?y) 


bjy 


(2.13) 


For  the  update  based  on  the  true  D  ,  yt  =  i>+2(x+-xc),  so  that 
(i+-ze)Ty(  =  (x+-xe)rZ>;2(x+-xe)  >  0  . 

For  the  update  based  on  the  barrier  function,  yb  =  )e,  and  we  have 

(■ x+-xc)Tyb  =  {x+-xc)TDc-1D-l(x+-xc)  >  0  . 

For  the  update  based  on  the  potential  function, 


y,  =  )r  - 


n(c  T xc-c  T x+)2 
(c  Tx+)2(cTxc) 


c  ,  so  that 
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(x+-xc)Ty  =  (x+- xc)T Dc  XD+  (x+-xc)  - 


n[cT{x+-xc)}3 
(c  Tx  +  )i{c  Txc) 


which  is  not  guaranteed  to  be  positive,  but  is  certainly  positive  if 
cTx+  <  xc  .  If  this  updating  strategy  is  being  used  and  the  above  is  not 
positive  at  some  step,  we  choose  to  use  the  update  based  on  the  true  D  ,  rather 
than  reject  the  update  altogether.  Thus  far  in  our  computational  testing,  the 
situation  in  which  ( x+-xc)Typ  <  0  has  not  arisen. 

The  updating  formula  given  by  equation  (2.10)  was  obtained  as  a  least- 
change  secant  update  for  approximating  D,  the  matrix  we  are  interested  in,  mak¬ 
ing  use  of  nonlinear  functions  whose  gradients  and  Hessians  involve  the  matrix 
Z)"1.  We  have  also  explored  the  use  of  the  updating  strategy  that  changes  the 
approximation  to  the  Hessian  of  the  nonlinear  function  as  little  as  possible  (rather 
than  changing  the  approximation  to  the  inverse  of  the  Hessian  as  little  as  possi¬ 
ble)  consistent  with  having  the  approximation  satisfy  the  secant  equation  (2.6)  or 
(2.7)  or  the  condition  (2.8).  Using  this  approach,  the  approximation  is  required  to 
satisfy 


D+D+T{x+-xc)  =  y  , 


and  analogous  to  (2.9),  we  express  this  requirement  as 

y  =  D?v 
v  =  D+  (x+  -  xc ) 

a.  _ i 

for  some  vector  v  G  1R”  .  Now  the  closest  matrix  Z>+ 
part  of  (2.14)  is 


(2.14) 


to  Dc  1  satisfying  the  first 


=  d;x  + 


(■ y-bcxv)v' 


T 

V  V 


(2.15) 


and  from  the  second  part  of  (2.14)  we  obtain 
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v 


_ yT(x+-*c) _ 

( Arr(x+ -  xt  ))t{dc~t(x+  -  xc )) 


* 

d;t(x+-xc ) . 


(2.18) 


Using  the  Sherman-Morrison-Woodbury  formula,  (2.15)  is  equivalent  to 

iv  ~Dc  y)(x+~ xc)T 


D+  =  Dc  + 


yr(*+-*c) 


(2.17) 


Use  of  the  update  given  by  (2.17)  requires  more  computation  and  more  storage 
than  is  required  if  the  update  given  by  (2.10)  is  used,  as  we  shall  discuss  in  sec¬ 
tion  3.  Furthermore,  numerical  experience  thus  far  has  indicated  better  perfor¬ 
mance  is  achieved  using  the  update  given  by  (2.10).  The  modified  algorithm  will 
be  referred  to  as  the  DMT-Karmarkar  algorithm. 


3.  Computational  issues  in  the  DMT-Karmarkar  algorithm 
3.1  Update  vectors 


The  computation  of  cp  requires  the  solution  of  a  linear  system  with 

A  A  *T» 

coefficient  matrix  BB  .  Since 


BBT  = 


M  a 


where  M  =  (AD)(AD)T ,  a=ADDTD~1e,  and  <7  =  e TD~1DDTD~1e,  we 
need  only  be  concerned  with  solving  a  linear  system  with  coefficient  matrix 
(AD)(AD)T .  We  first  discuss  the  computation  of  cp  using  the  update  to  D  given 
by  (2.10)  and  (2.13).  Referring  to  the  updating  formula  (2.10), 


AD+  =  ADC  - 


{ADcv)v‘ 

V  T  V 


=  ADr 


I  - 


w 

T 

V  V 


Hence, 
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(AD^AD+)t  =  AD,\l-  ^-\ (AD, )’ 

l  V1  v) 


=  (Ab')(Af>')T  - 


(ADcv)(ADev)T 


Now  letting  Mc  and  M+  denote  (ADC)(ADC )T  and  ( AD+)(AD+)T ,  respectively, 

A  m 

and  setting  w  —  ADcv  and  0  —  v  v,  the  Sherman-Morrison- Woodbury  formula 


gives 


M =  \MC  -  SL 


,  M~1wivTMcl 

=  M-1  +  — ^ - p- 

0-wT  Mc  lv) 


Setting  t  =  Mc  lw  and  ^  —  0  -  wTMc  lw  ,  this  is 


In  summary,  at  each  step  of  the  algorithm  we  save  four  vectors  and  two 
scalars  :  v,u  £  Rn,  t u,t  £  IRm,  and  0 ,7  £  1R  as  follows  : 


=  (x+-x*)Ty 

A 

«  =  x+  -  xc  -  Dc  v  , 


Dc  y  ,  w  =  ADcv  , 


(  -  \(ADc)(Abc)T\-'w  , 


(3.1.1) 


0  =  vTv  , 

The  computation  of  c.  is  carried  out  using 


=  0  -  wT t 


Dk=D0+  -i-i-  + 
P\ 


(3.1.2) 


=  /  + 


T 

4 


T 

1  +  -“H  [(^o)(AZ>o)rrS(3-1.3) 


so  that  the  only  factorization  needed  for  £  steps  of  the  algorithm  is  of  the  initial 
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matrix  ( AD0)(AD0)T  .  Notice  also  that  vector  or  pipeline  architectures  are 
immediately  applicable  to  the  required  computations. 

We  are  indebted  to  M.J.  Todd  (private  communication)  for  pointing  out  to 
us  that  computational  savings  may  be  realized  if  the  updates  to  A/-1  are  saved  in 
factored  form.  To  this  end,  suppose  that  (ADC\ADC)T  =  LcLj.  Then 
l(ADc)(ADc)T}-'  =  L~TL~l  m  NCN ?  ,  and 


M 


-i  _  „  ■  NcNcT™TNcNeT 


NCNC 1  + 


0-  wtNcNctw 


Setting  q  =  N?w  and  7  —  0  -  qTq  +  [0(0  -  qTq)f‘ ,  we  obtain 

M;1  =NcNct  + 


NcqqTN J 


0-  qTq 


so  that  the  update  to  the  Cholesky  factor  is  given  by 

,T 


N,=N.\I  + 


Using  this  approach,  one  would  save,  at  each  step  of  the  algorithm,  three 
vectors  and  two  scalars  :  v,u  £  ]Rn,  q  £  IRm,  and  /?,7  £1R  as  follows  : 

Djy  ,  q  =  NjADc  v  , 


(x+-xc)Ty 


(Djy)T(£>J  y) 


u  =  -  xc  -  Dcv  , 


0=  vTv  , 


7  =  0-  qTq  +  [0(0-  • 


IS 


The  computation  of  cf  would  be  carried  out  using  (3.1.2)  and 

[(AD„)(AD,,)T}-' =  NtN?  (3.1.4) 

where 


Nt  =  Nr 


I  + 


9i9i 


h 


/  + 


T  \ 
9292 
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I  + 


T 

Qk  9* 

— T  1 

Ik 


Saving  the  updates  in  factored  form  would  save  storage  of  w  ,  and  the  computa¬ 
tion  of  q  rather  than  t  would  save  n  floating  point  operations,  where  n  is  the 
number  of  nonzeros  in  the  Cholesky  factor  of  ( AD0)(AD0)T  .  The  computation 
(3.1.4),  however,  would  require  mk  more  floating  point  operations  than  (3.1.3)  at 
the  Ar th  step,  so  that  the  amount  of  computation  saved  would  vary,  depending  on 
sparsity.  The  implementation  for  which  we  give  numerical  results  in  section  4 
uses  updates  (3.1.1). 

If  the  updating  formula  given  by  (2.17)  supported  by  (2.16)  is  used,  the  com¬ 
putation  of  Cp  requires  twice  as  much  storage  and  nearly  twice  as  much  work  as 
this  computation  using  updating  formula  (2.10)  supported  by  (2.13).  This  is 
because  the  use  of  (2.17)  leads  to  a  rank-two  update  to  obtain  {AD +\AD +)T 
from  ( ADc){ADc)t  at  each  step.  At  each  step  of  the  modified  algorithm  using 
the  update  (2.17),  the  eight  vectors  s,v,u,t  G  R”  and  p,q,d,g  E  Rm  and  the 
four  scalars  /Vy,7/,  and  p  are  saved  : 


s  =  x+-xc, 


v 


T 

y  * 


{d;ts)t{d;ts) 


* 

Dc~Ts , 


u  =  v  - Dcy, 

t  =  b;lu, 


P  =  Au, 
q  =  ADC  s , 

d  =  M-\q  +  ^p), 
9  =  A  Clte  P  ’ 
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0  =  VTs, 

1  =  vTv, 


where  Mc  =  ( ADC)(ADC)T  , 


V  =  P  +  PTd, 

p  =  0  +  qTg, 


=  /- 


9k-  \9k-\ 
Pk-i 


I  - 


9i9i 

Pi 


dk-\Vk-\ 


Vk-i 

1 
i 


d\PT 


*h 


{(AD0)(AD0)TY 


and 


mk-yt 


=  \I- 


dkPk 


Vk 


I  - 


I  - 


I  - 


9i9j 

Pi 


9k-\9k-\ 
Pk-l 
d  vT 


dk-lP?-l 

Vk-l 


The  approximation  Dk  is  given  by 


At  =  + 


«i*i 

01 


+ 


+ 


w 


0k  ' 

The  computation  of  v  requires  D~Ts.  The  matrix  Af1  is 

hvk 


A r 1  =  ^o1 


livi 


ik 


Using  updating  formula  (2.17),  it  is  still  the  case  that  the  only  matrix  factoriza¬ 
tion  needed  for  k  steps  of  the  algorithm  is  of  the  initial  matrix  (AD0)(AD0)T  . 


3.2  Restarting  strategy  for  the  DMT-Karmarkar  algorithm 

In  practice,  there  is  a  limit  on  the  number  of  update  vectors  one  is  willing  to 
store.  We  have  implemented  two  ways  of  proceeding  when  the  allotted  storage 
has  been  used.  The  option  we  recommend  is  to  restart,  treating  the  current  point 
as  the  initial  point.  This,  of  course,  requires  another  matrix  factorization.  In  our 
implementation,  the  restart  strategy  includes  restarting  when  the  approximation 
appears  not  to  give  a  direction  of  sufficient  decrease  as  well  as  when  the 
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maximum  number  of  update  vectors  has  been  saved.  We  also  plan  to  try  a  res¬ 
tart  strategy  based  on  the  partial  refactorization  suggested  by  Karmarkar  [1984] 
and  studied  by  Shanno  [1985b].  A  second  option  that  we  have  considered  but  do 
not  recommend  is  to  discard  all  previous  updates,  and  treat  the  current  point  as 
the  first  point;  that  is,  to  replace  the  collection  of  updates  with  a  single  update  as 
if  the  current  point  had  been  reached  in  one  step  from  the  initial  point.  One 
may  also  combine  the  two  strategies,  electing  to  perform  an  additional  matrix 
factorization  only  when  it  is  determined  that  progress  is  not  being  made.  Using 
the  pure  strategy  of  discarding  updates  means  that  only  the  initial  matrix  is  fac¬ 
tored.  This  can  yield  the  solution  to  the  problem  at  modest  expense  and  using 
very  limited  storage,  but  convergence  in  this  case  is  not  guaranteed.  We  recom¬ 
mend  the  restarting  strategy  over  the  discarding  strategy,  both  for  its  sound 
theoretical  basis  and  for  its  superior  performance  in  practice.  Use  of  the  restart¬ 
ing  strategy  results  in  an  algorithm  that  retains  the  polynomial  worst  case  time 
bound  of  the  Karmarkar  algorithm,  since  our  step  acceptance  criterion,  to  be  dis¬ 
cussed  below,  requires  reduction  of  the  potential  function  by  at  least  the  constant 
amount  at  each  step  that  is  guaranteed  for  the  Karmarkar  algorithm,  and  clearly 
one  can  always  get  this  reduction  by  restarting. 

A  natural  question  to  ask  is  whether  it  is  advantageous  to  use  update  vec¬ 
tors  at  all  in  the  DMT-Karmarkar  algorithm.  A  reasonable  alternative  strategy 
within  the  framework  we  have  established  in  (2.1)  is,  as  mentioned  by  Gill,  Mur¬ 
ray,  Saunders,  and  Wright  [1986],  to  retain  the  same  approximation  to  D  for 
several  iterations  before  recomputing  a  factorization.  In  an  algorithm  with 
periodic  restarts,  one  would  have  /),  =  D0  for  0  <  *  <  Ar— 1  for  some  k,  and 
then  treat  the  Ar th  iterate  as  the  initial  point.  We  have  also  explored  the  use  of 
this  strategy. 
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S.S  Linesearch  and  step  acceptance  criterion 

Our  implementation  of  the  DMT-Karmarkar  algorithm  includes  a  linesearch 
with  a  three  faceted  step  acceptance  criterion  as  a  way  of  selecting  the 
steplength.  The  step  acceptance  criterion  ensures  that  our  algorithm  with  res¬ 
tarts  retains  the  polynomial  time  bound  of  the  unmodified  Karmarkar  algorithm 
but,  rather  than  restricting  steps  to  have  a  predetermined  fixed  length,  allows 
taking  longer  steps  when  it  appears  advantageous  to  do  so.  We  interpret  failure 
to  find  an  acceptable  step  after  a  specified  number  of  trial  steps  as  an  indication 
that  a  restart  is  needed.  Our  first  trial  steplength  is  .99  of  the  distance  to  the 
edge  of  the  simplex;  our  third  trial  steplength  is  .99  of  the  radius  of  the  largest 
sphere  that  can  be  inscribed  in  the  simplex,  and  our  second  trial  steplength  is 
midway  between  the  first  and  third.  Beginning  with  the  fourth  trial  step,  each  is 
half  as  long  as  the  previous  trial  step.  Our  step  acceptance  criterion  is  a  combi¬ 
nation  of  the  Karmarkar  criterion,  the  Goldstein-Armijo  condition,  and  reduction 
of  the  linear  objective  function.  The  Karmarkar  criterion,  which  is  used  to  prove 
convergence  of  the  algorithm  in  polynomial  time,  requires  at  least  a  constant 
amount  of  reduction  in  the  potential  function  (2.3)  at  every  step: 

/(«-«*,)</(.)-*  (3.3.1) 

=  n  log  c  Txc  -  6, 

where  6  is  the  minimum  reduction  in  the  potential  function  that  is  guaranteed  for 
each  step  of  the  Karmarkar  algorithm.  The  Goldstein-Armijo  condition  (see,  for 
example,  Dennis  and  Schnabel  [1983]),  requires  that  the  average  rate  of  decrease 
in  the  potential  function  from  x  c  to  x  +  be  at  least  a  prescribed  fraction  (X)  of 
the  initial  rate  of  decrease  in  that  direction: 

/  (*  +)  <  /  (*  c)  +  x  V/  {i  c)T(i  +-i  c)  ; 

that  is, 
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f  (c-ac  )  <  n  log  cTxe  -  a 


Xn 


T 

c1  x. 


{Dc)Tc. 


(3.3.2) 


In  order  to  determine  whether  the  trial  step  would  result  in  a  reduction  of  the 
linear  function  cT x  ,  we  consider  the  objective  value  of  the  image  of  the  trial 
step  in  the  untransformed  space,  and  determine  whether  the  trial  step  satisfies 


T 

*  nr  - — 


cx 


n 


c  tDx  + 


■c T  Dx  ,  <  c 


(3.3.3) 


The  lower  bound  6  on  reduction  in  the  potential  function  obtained  by  Karmarkar 
is  dependent  on  the  size  of  the  problem  as  well  as  the  steplength,  and  ensures 
reduction  of  the  potential  function  at  each  step  by  at  least  .1  for  a  step  in  the 
direction  -cp  of  length  equal  to  one-fourth  the  radius  of  the  largest  sphere  that 
can  be  inscribed  in  the  simplex,  provided  the  number  of  variables  in  the  problem 
is  at  least  21.  Since  our  test  problems  are  larger  than  that,  we  use  6  =  .1  . 


Our  real  objective  is  minimization  of  the  linear  function  c  Tx  .  While  the 
linear  function  is  not  transformed  to  a  linear  function  under  projective  transfor¬ 
mation,  straight  lines  do  map  to  straight  lines  under  projective  transformation,  so 
that  the  level  sets  of  the  linear  function  are  mapped  to  flats  in  the  transformed 
space.  Therefore,  if  the  direction  we  are  considering  gives  descent  on  the  linear 
function,  it  seems  sensible  to  take  the  longest  possible  step  in  that  direction  that 
satisfies  (3.3.1).  In  the  event  our  direction  is  not  a  descent  direction  for  the  linear 
function,  we  guard  against  taking  a  step  too  long  relative  to  the  amount  of 
reduction  achieved  in  the  potential  function  by  requiring  our  step  to  satisfy 

(3.3.2)  as  well  as  (3.3.1).  That  is  to  say,  we  accept  the  step  as  soon  as  the  trial 
step  satisfies  (3.3.1)  and  either  (3.3.2)  or  (3.3.3).  For  an  iteration  at  which  a  res¬ 
tart  has  been  done,  if  a  longer  step  does  not  satisfy  (3.3.1)  and  either  (3.3.2)  or 

(3.3.3) ,  we  accept  a  step  in  the  direction  -cp  of  length  equal  to  one-fourth  the 
radius  of  the  largest  inscribed  sphere.  It  will  always  satisfy  (3.3.1). 
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4.  Numerical  results 

Preliminary  testing  using  an  experimental  implementation  of  the  DMT- 
Karmarkar  algorithm  has  given  very  encouraging  computational  results.  We 
have  used  eight  test  problems.  The  first  is  a  problem  with  one  hundred  seven 
variables  and  sixty-seven  constraints  that  we  obtained  from  Shell  Development 
Company;  the  second  through  eighth  are  test  problems  from  the  Systems  Optimi¬ 
zation  Laboratory  at  Stanford  University  that  we  obtained  through  netlib  (see 
Gay  [1985b]).  The  number  of  variables  shown  for  each  problem  includes  slack 
variables  introduced  to  transform  inequality  constraints  to  equalities,  but  does 
not  include  the  two  additional  slack  variables  discussed  in  Appendix  I.  Results 
obtained  solving  these  problems  using  various  options  in  the  DMT-Karmarkar 
algorithm  are  summarized  in  tables  4.1  through  4.8.  For  each  problem,  a  feasible 
starting  point  x0  was  obtained  using  the  technique  advocated  by  Karmarkar, 
which  we  include  in  Appendix  I,  from  the  initial  point  (1,  •  ••  ,1  )T  .  All  options 
were  then  run  from  the  same  feasible  starting  point.  The  number  of  steps  to  con¬ 
vergence  shown  in  the  tables  does  not  include  steps  required  to  obtain  the  start¬ 
ing  point.  The  first  option  shown  for  each  problem  is  restarting  every  step, 
which  means  that  the  Karmarkar  algorithm,  modified  only  by  the  incorporation 
of  a  linesearch,  is  used.  In  each  instance,  restart  after  k  updates  really  means 
that  a  restart  is  done  after  saving  at  most  k  updates;  that  is,  at  least  as  often 
as  every  k  + 1  steps.  A  new  factorization  is  done  earlier  whenever  the  linesearch 
is  unsuccessful.  Similarly,  restart  after  k  steps  (using  no  updates,  but  retaining 
the  same  approximation  to  D  between  factorizations)  means  restart  after  at 
most  k  steps.  To  test  the  DMT-Karmarkar  algorithm  against  the  Karmarkar 
algorithm,  we  solved  the  test  problems  to  only  limited  accuracy.  We  expect  that 
relative  times  needed  to  compute  more  accurate  solutions  would  be  similar  to 
those  given  in  tables  4.1  through  4.8.  For  all  problems,  the  stopping  criterion 
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was  c  ^ x  <  lO^c  T Xq  .  This  gave  three  to  four  digits  of  accuracy  in  the  optimal 
objective  values  for  all  of  the  test  problems  except  ISRAEL,  where  only  one  digit 
of  accuracy  in  the  optimal  objective  value  was  obtained.  Subsequently,  we  con¬ 
tinued  the  solution  of  ISRAEL.  We  obtained  four  digits  of  accuracy  in  the 
optimal  objective  value  at  the  cost  of  twelve  additional  steps  and  factorizations 
using  the  Karmarkar  algorithm,  seventeen  steps  with  nine  factorizations  restart¬ 
ing  after  one  update,  twenty-eight  steps  with  seven  factorizations  restarting  after 
three  updates,  and  thirty-four  steps  with  six  factorizations  restarting  after  five 
updates. 

In  tables  4.1  through  4.8  the  update  used  is  identified  as  follows  : 


update  number  updating  formula  righthand  side  of 

secant  equation 

1 

2.10' 

yb  of  (2.61 

2 

2.17 

yb  of  (2.6) 

3 

2.17 

yp  of  2.7 
yp  of  2.7) 

4 

2.10 

5 

2.10 

yt  of  (2.8) 

6 

2.17 

Vt  of  (2.8) 

For  options  involving  the  use  of  updates,  we  show  only  the  results  obtained 
using  the  most  successful  updating  formula.  In  general,  better  performance  was 
achieved  using  updating  formulas  1,  4  and  5,  the  least  change  updates  to  D 
based  on  the  barrier  function,  the  potential  function,  and  on  the  true  D,  respec¬ 
tively,  than  was  achieved  using  updating  formulas  2,  3  and  6,  the  updates  that 
result  from  the  strategy  of  least  change  to  the  Hessians  (t.e.  least  change  to  jD_1) 
of  the  functions  used  to  obtain  the  approximations.  The  most  successful  updat¬ 
ing  formula  was  the  simplest,  the  least  change  update  to  D  based  on  the  true  D, 
identified  as  update  5.  In  cases  where  two  or  more  updating  formulas  produced 
identical  results  in  number  of  steps  and  factorizations  needed  to  obtain  the  solu- 
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tions,  times  given  are  for  updating  formula  5. 

We  have  obtained  timings  solving  the  test  problems  on  a  Pyramid  90x,  using 
the  Unix  operating  system  OSx  version  2.5.  Times  are  the  sum  of  CPU  times 
attributed  to  the  user  and  to  the  operating  system  in  seconds,  rounded  to  the 
nearest  second.  The  last  column  of  each  table  contains  normalized  times,  the 
ratio  of  the  time  required  to  obtain  the  solution  using  the  given  option  to  the 
time  required  using  the  Karmarkar  algorithm  (with  linesearch). 

For  each  of  the  test  problems,  use  of  the  DMT-Karmarkar  algorithm  with 
periodic  restarts,  and  using  rank-one  updates  to  approximate  D  at  each  step, 
results  in  obtaining  the  solution  with  (generally  but  not  monotonically)  progres¬ 
sively  fewer  matrix  factorizations  as  the  number  of  updates  allowed  between  fac¬ 
torizations  increases.  Of  course,  the  number  of  iterations  required  to  obtain  the 
solution  increases  as  the  number  of  factorizations  decreases.  Timings  on  the  test 
problems  indicate  an  overall  reduction  in  the  amount  of  work  can  be  achieved 
using  the  DMT-Karmarkar  algorithm,  compared  to  the  Karmarkar  algorithm,  for 
all  problems  except  our  smallest  test  problem,  AFIRO  (Table  4.2).  Timings  are 
accurate  only  to  within  about  one  second,  so  that  little  significance  should  be 
given  to  relative  times  for  solving  this  very  small  problem  that  takes  only  two  to 
three  seconds.  The  largest  percentages  of  savings  in  computational  effort  appear, 
as  one  would  expect,  in  the  larger  problems,  ISRAEL  (Table  4.7)  and  BRANDY 
(Table  4.8). 

The  strategy  of  periodic  restarts  with  no  updates  used  may  be  advantageous 
if  only  a  few  steps  are  taken  between  restarts.  The  increase  in  number  of  steps 
required  is  generally  greater,  and  the  number  of  factorizations  saved  is  generally 
smaller,  than  occurs  when  update  vectors  are  used,  but  the  computation  of  steps 
that  require  neither  factorizations  nor  updates  is  very  inexpensive. 
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Test  problem:  SHELL:  107  variables,  67  constraints 

Option 

Update 

Steps 

F  actorizations 

Time 

Time/Time(K) 

Restart 

8 

8 

17 

1.00 

5 

12 

6 

16 

.94 

5 

14 

5 

15 

.88 

K53HE5BSsB 

5 

18 

4 

16 

.94 

after  9  updates 

5 

28 

3 

21 

1.24 

after  1  step 

mam 

13 

7 

17 

1.00 

1 

ii 

19 

7 

19 

1.12 

IlSSEES^Bi 

IMM 

21 

7 

20 

1.18 

Table  4.1 


Test  problem:  AFIRO:  51  variables,  27  constraints 


Option 

Update 

Steps 

F  actorizations 

Time 

Time/Time(K) 

Restart 

^EE53E39HHI 

7 

7 

2 

1.00 

B3BBSS9II 

IB591 

11 

6 

3 

1.50 

5 

12 

5 

3 

1.50 

after  3  updates 

5 

16 

4 

3 

1.50 

after  7  updates 

5 

24 

3 

5 

2.50 

after  1  step 

11 

6 

3 

1.50 

mmm 

16 

6 

3 

1.50 

IKESKESSi 

IBM 

21 

6 

4 

2.00 

Table  4.2 


Test  problem:  SHARE2B:  162  variables,  96 

const  rai 

nts 

Option 

Update 

Steps 

Factorizations 

Time 

Time/Time(K) 

Restart 

every  step 

9 

9 

26 

1.00 

5  ... 

14 

7 

24 

.92 

after  3  updates 

5 

22 

6 

28 

1.08 

after  6  updates 

33 

5 

36 

1.38 

after  1  step 

mmm 

16  ™ 

8 

27 

1.04 

mam 

22 

8 

29 

1.12 

IMISBMMSSM 

warn 

41 

7 

35 

1.35 

Table  4.4 
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Test  pi 

■oblem:  SHARElB:  253  variables,  117  constraints 

Option 

Update 

Steps 

F  actorizations 

Time 

Time/Time(K) 

Restart 

19 

19 

248 

1.00 

5 

25 

13 

198 

.80 

5 

31 

11 

175 

.71 

5 

34 

9 

158 

.64 

5 

39 

8 

155 

.63 

after  5  updates 

5 

61 

11 

222 

.90 

5 

49 

7 

163 

.66 

after  0  updates 

5 

58 

_  6 

166 

.67 

after  1  step 

mmm 

28 

14 

193 

.78 

IMHI 

37 

13 

192 

.77 

\mmm 

44 

12 

190 

.77 

Table  4.5 


Test  problem:  BEACONFD:  295  variables,  173  constraints 


Update  Steps  Factorizations  Time  Time/Time(K) 


273 

1.00 

227 

.83 

215 

.79 

200 

.73 

212 

00 

• 

227 

.83 

222 

.81 

232 

.85 

241 

.88 

Test  pi 

•oblem:  BR 

ANDY:  292  variables,  182  constraints 

Option 

Update 

Steps 

Factorizations 

Time 

Time/Time(K) 

Restart 

mmm 

12 

12 

276 

1.00 

wnm 

17 

9 

223 

.81 

lESSHSJSrai 

5 

20 

7 

189 

.68 

5 

24 

6 

176 

.64 

5 

32 

5 

175 

.64 

after  10  updates 

5 

41 

4 

ms 

.65 

after  1  step 

!■■ 

20 

10 

244 

.88 

iX2&9E£5&9H 

imm 

25 

9 

229 

.83 

■ 

33 

9 

240 

.87 

Table  4.8 


Appendix  I 


The  standard  linear  programming  problem 

.  .  .  —T— 

minimize  c  x 

subject  to  Ax  —  b 

x  >  0, 

where  c,x  6  TRn~2,  b  6  IR"-1,  and  A  6  r(*-»)x(»-2)  can  be  transformed  into  the 
linear  programming  problem 


minimize  cTx 

subject  to  Ax  —  0 
T 

e  x  =  n 
x  >  0, 


(A.1) 


where  c,x,e  e  IR",  A  €  IRmXB,  and  e=(l,  •  •  •  ,1)T  as  follows.  First,  add  a 
slack  variable  x*^  ,  express  the  requirement  Ax  —  b  as  Ax  —  ib_16  ,  and  add 
a  constraint  requiring  =  1  .  The  equality  constraints  can  now  be  stated  as 


A 


*1 

o' 

*n-  2 

..  - 

0 

1 

1 _ 

1 

Now  assume  we  have  a  bound  B  on  the  sum  of  the  variables  : 

n- 1 

<B 

t=l 

and  introduce  a  second  slack  variable  so  that 


E  *;  =  ■£> 

i  =  i 

Using  this  condition,  the  constraint  equation  we  have  added  can  be  written  as 


27 


*-■  =  J  t* 

Hence,  the  equality  constraints  have  been  put  into  homogeneous  form: 


Now  scale  the  variables  so  that  their  sum  is  n  :  x  —  —  x  .  Since  the  number 

B 

of  variables  has  been  increased  by  two,  two  zero  components  are  appended  to  c 
to  obtain  c  . 

If  the  optimal  value  of  the  objective  function  is  not  zero,  but  is  known  to  be 
f*  ,  the  objective  function  can  easily  be  transformed  to  an  objective  function 
whose  minimum  value  is  zero.  Minimizing  cTx  is  equivalent  to  minimizing 
c  x  -  /  *  ,  and  this  second  objective  function  is  equivalent  to 


so  that 


E  *«■  f 


i=l 


For  our  test  problems,  the  optimal  value  of  the  objective  function  was 
known,  so  that  only  the  shift  in  c  as  above  was  necessary.  If  the  minimum  value 
of  the  objective  function  is  not  known,  Karmarkar  [1984]  suggests  the  use  of 
what  he  calls  a  sliding  objective  function.  This  involves  attempting  to  solve  the 

problem  with  objective  function  c  T x  with  c  =  c  -  —  e  ,  where  /  is  an  estimate 

n 

of  f  *  that  is  refined  as  the  solution  progresses.  Todd  and  Burrell  [1985]  have 
suggested  using  the  dual  of  problem  (A.l)  to  obtain  an  estimate  /  that  is  a  lower 


S8 


bound  on  /  *  and  that  is  refined  at  each  iteration. 

As  noted  by  Karmarkar  [1984],  given  any  initial  point  x  G  Rn  ,  x  >  0  ,  and 
assuming  the  linear  programming  problem  is  feasible,  an  initial  feasible  point  for 
the  linear  programming  problem  (A.l)  can  be  obtained  by  solving  (for  x  G  Rn+1) 
the  problem 

minimize  xn+1 

subject  to  Ax  =  0 

c  T  x  =  n  +1 
x  >  0, 

where  A  G  Rm  x(n+1)  has  as  its  first  n  columns  the  matrix  A  and  as  its  last 
column  the  vector  -Ax  .  A  feasible  starting  point  for  problem  (A.3)  is 
(*i,  !)r  • 
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